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Abstract. This paper deals with a problem from discrete-time robust 
control which requires the solution of constraints over the reals that con- 
tain both universal and existential quantifiers. For solving this problem 
we formulate it as a program in a (fictitious) constraint logic program- 
ming language with explicit quantifier notation. This allows us to clarify 
the special structure of the problem, and to extend an algorithm for 
computing approximate solution sets of first-order constraints over the 
reals to exploit this structure. As a result we can deal with inputs that 
are clearly out of reach for current symbolic solvers. 



1 Introduction 

A discrete time system can be described by state-space equations of the form 
Xk+i = f{xk,Uk,'Wk), where k is the discrete time, / is a non-linear real func- 
tion, Xk is the state vector at time fc, Uk is the control vector which can be 
chosen arbitrarily in a set Uk, and Wk is the perturbation vector which cannot 
be influenced by us but remains inside a known set Wk- In this paper we deal 
with the problem of computing the set of state vectors Xq for which we can set 
the controls in such a way that future state vectors Xi, . . . ,Xn will belong to 
certain sets Xi, . . . , X„ chosen by us. 

This problem is closely related to the problem of characterization of viability 
sets involved when studying the evolution of macro-systems arising in biology, 
economics, cognitive sciences, games, and similar areas, as well as in nonlinear 
systems of control theory. A good reference is the book of Aubin |^ . 

When trying to solve this problem one immediately runs into constraints 
that contain a large number of universally and cxistentially quantified variables. 
In this paper we report on ongoing research on a method for solving such con- 
straints. This method is already able to compute (approximate) solutions that 
are far too hard to compute for current symbolic solvers (e.g., QEPCAD |p^). 

We proceed by giving a recursive formulation of the problem which can be 
read as a program in a constraint logic programming language with explicit 
quantifier notation. This formulation will allow us to clarify the special structure 
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of the resulting first-order constraint (i.e., formula in the first-order predicate 
language with predicate symbols = and <, function symbols + and x, and 
rational constants, all of them with their usual interpretation). We exploit this 
structure by breaking the constraint into parts, where each part corresponds to 
one stage of the system, and these parts are glued together by a syntactic entity 
called "function inversion quantifier" . We show how to extend the method of 
approximate quantified constraint solving (AQCS |32| , |30|l ), as developed by one 
of the authors, to deal with such function inversion quantifiers. This will allow 
us to solve non-trivial instances of the mentioned problem from discrete-time 
robust control. 

The resulting approach is also applicable to several other problems from 
discrete-time control, and even to the general case of first-order constraints 
with composition structure, that is, first-order constraints that have the form 
C{fn{- ■ ■ fi{x) . . .)). Since defining objects in hierarchies is ubiquitous in the 
human modeling and problem-solving process, and since such a composition 
structure very often is (implicitly) encoded in programs in CLP languages, we 
believe that according support in constraint solvers will become more and more 
important. 



2 Problem Definition 

A discrete-time system with control and perturbation (or short a system) is a 
tuple (/, n, X, U, W), where 

- / : IR"+"+«' ^ H"^, 

- n e IN, 

- AC {l,...,n} xWi", 

- U C {l,...,n} X WC\ and 

- W C{l,...,n} X R"'. 

The function / is called transition function, the predicate A allowed state, the 
predicate U allowed control, the predicate W possible perturbation. Intuitively 
(see Figure |I|), such a system models a process which starts with some value 
Xq (the initial state) and then applies the function / n-times to Xq, where at 
each application additional user-definable input (control) s.t. U{k,Uk), and 
perturbation Wk s.t. W{k, Wk) influence the function / — resulting in the state at 
stage k. In this paper we solve the problem of computing the robust feasible initial 
set: Given a discrete-time system with control and perturbation find an initial 
state such that for all following stages the state is allowed (i.e., the predicate A 
holds). 

Given a system S = (/, n. A, U, W) we can formalize this using the following 
predicate, which models the states x at stage k for which we can apply input 
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Fig. 1. Discrete Time System 



such that for all possible perturbations the state at all future stages is allowed: 



Cs{x,k) :< — >A{x,k)A 

k < n ^ 

3u U{u,k) A 
Vw W{w,k) 



(1) 



Cs{fix,u,w),k + 1) 



Our problem is to compute the robust feasible initial set Cs{xo,0). After 
expanding this query and writing the predicates U and W as sets, we get: 



Ao{xo) A 3uo G Uo e Wa 

Ai{f{xo,UQ,wo)) A 3ui £ Ui ywi e Wi 

A2{f{f{xo,uo,wo),Ui,Wi)) A 3u2 e U2 e W2 

An{f{f{- . . f{Xo,Uo,Wo),Ui,Wi),U2,W2)) 



(2) 



This is an expression that is definitely too hard to solve for current symbolic 
solvers like QEPCAD (for which the border between solvable and unsolvable 
problems is around 3-5 variables), and in its raw form also for the approximate 
solver developed by one of the authors ^2 30 . Thus we have to exploit the 
special structure of such a constraint in order to be able to solve it. 



3 Main Idea 

Observe that in the recursive definition of in Expression we only have a 
fixed number of variables in each recursion step. When expanding the definition, 
each step maps the newly introduced variables into a variable vector of lower 
dimension by applying the function /. In the other direction, when collecting 
results, we can first compute the solution set of step n, plug the result into step 
n — 1, compute this solution set, and so on. 

Here we have to solve a first-order constraint in a + u + w variables for each 
step. But even if we could compute the solution set of step n by a symbolic 
solver, the size of the result will blow up when plugging it into the preceding 
steps n — 1, n — 2, and so on. The reason is that the complexity of real quantifier 
elimination results from the size of its quantifier-free output [ ^Jl3| . So we use 
the same idea to compute approximate solutions instead. 



For this we abstract from our concrete example and assume a constraint C for 
which we want to compute the solution set of C{f{x)), where the dimension of 
X is higher than the dimension of f{x). However, we want to avoid substituting 
f{x) into the constraint because that would blow up the complexity by forcing 
us to solve a constraint with the additional variables x. 

Thus we try to deal with such a situation in a special, more efficient way. 
For this we write it by using Junction inversion quantifiers^ that is, syntactical 
entities of the form [f{x) = z] C{z). Semantically such a constraint is equivalent 
to C{f{x)). Operationally it will be handled in a different, more efficient way. 
Now our problem reads as follows: 

Aq{xq) a 3uo e Uq Vwo e Wo [f{xo, uo,wo) = xi] 
Ai{xi) A 3ui e Ui Vwi G Wi [f ixi,ui,wi) ^ X2] 

A2{X2) A 3U2 e U2 yW2 e W2 [f ix2, U2,W2) = X3] (3) 
An {Xji^ 

Observe that in this constraint each sub-constraint beginning at a certain 
line has exactly one free- variable vector Xi. Now we can view our problem of 
computing the robust feasible initial set as the problem of solving n constraints 
that are glued together by function inversion quantifiers. For 1 < i < n, the j-th 
of these constraints only contains the variable vectors Xi, Ui, and Wi. 

For solving these constraints we use the method for approximate quantified 
constraint solving (AQCS) as developed by the first author In the rest 

of the paper we will describe how to glue together several instances of this solver 
by implementing function inversion quantifiers. However, the scope of this paper 
will only allow us to explain these details of AQCS that are absolutely necessary 
for understanding how we deal with function inversion quantifiers. 

4 Function Inversion in Quantified Constraint Solving 

In approximate quantified constraint solving we approximate solution sets by 
functions that assign "true" to elements that are guaranteed to be in the solution 
set, "false" to elements that are guaranteed to be out of the solution set, and 
that assign "unknown" otherwise: 

Definition 1. An approximate set on i? C IR" is a function in B {T, F, U}. 

The error err(5) of an approximate set S on B is the volume of the x such that 
S{x) ~ U. An approximate set S on B is an approximation of a set S ^ B iff 
for all X G B, S{x) = T implies that x E S, and S{x) = F implies x ^ S. 

We delay the description of a concrete computer representation of approximate 
sets to the end of this section. 

The basic idea of the algorithm for approximate quantified constraint solv- 
ing is to define for every syntactic element p (e.g.. A, V) a propagation 
function prop^^ that computes an approximation of the solution set of a con- 
straint of the form p{Ci , . . . ,Cn) from approximations of the solution set of each 



sub-constraint Ci, . . . , C„ |3l|. For example, for the constraint 3x + y"^ < 1 
the propagation function propg takes an approximation of the solution set of 
x'^ + < 1 and returns an approximation of the solution set of 3x x'^ + y'^ < 1. 
We follow this approach by simply defining such a propagation function also 
for function inversion quantifiers. This means that we show how to compute an 
approximation of the solution set of [f{x) — z]C{z) (i.e., C{f{x))) from an 
approximation of the solution set of C{z). In a first attempt we define for an 
approximate set S on B, prop[^(^-,^j,] (5) to be Xx e B.S{f{x)'j^. We will see 
later how to implement this in detail. 

Now a naive approach would apply the idea from Section ^ by computing an 
approximate solution set of stage n, then using this information to compute an 
approximate solution set of stage n — 1, and so on (this corresponds to a compu- 
tation that follows the arrows in reverse order in Figure |l|). The problem is that 
here we do not know how small we have to make the error of these approximate 
solution sets in order to reach a certain desired output error. Furthermore we do 
not know which parts of these approximate solution sets are needed for the end 
result. Thus we compute all of them on demand. 

This is also how approximate solution sets are computed within AQCS. For 
this it uses a function refine that takes a constraint C and a set i?, and returns 
an approximation of the solution set of C on i? (an exact discussion of how 
large the error of this approximation is allowed to be is beyond the scope of 
this paper). The argument B is chosen on demand within the algorithm — it 
is usually just a small part of the solution set we want to compute, and during 
the algorithm execution the limit of its volume goes to zero. So we have to 
extend this refinement function for the case when the outermost symbol of C is 
a function inversion quantifier: 

refinc([/(a;) z]C,B) = 

prop[/(^)=^](refine(C, f{B))) 

The resulting overall function is recursive, the base case is reached when 
called with an atomic constraint (i.e., an equality or inequality). We will deal 
with the problem of how to compute f{B) (i.e., the range of / on B) later. 

However, this solution does not yet suffice to solve the complexity problem. 
The reason is that for each input set B it computes an approximate solution 
set of C on f{B). In other words, it substitutes f{x) into C at the level of 
solving instead of at the syntactic level, and we do not gain anything in terms 
of complexity. 

This problem comes up because we never represent the approximate solution 
set of the sub-constraint C of [f{x) = z]C explicitly. This means that, even if 
for many different x the value of f{x) is the same, we recompute the solution set 
of C for different x that result in the same f{x) by new calls to the refinement 
function. We solve this problem by remembering already computed parts of the 



^ The notation Xx G B.f{x) denotes a function that takes an argument x £ B and 
returns f{x) |5[. 



solution set of C. We do this by wrapping the above caU to the refinement func- 
tion into the following memoization function which keeps the current knowledge 
about the solution set of C by caching it in the global variable S, which we 
initialize by the everywhere unknown approximate solution set Aa;GlR.U: 

refineMemo(C, B) = 

if err(restr(S', B)) > e then 

S ^ embed(S', refine(C, choose(S', B))) 
return restr(S', B) 

Here e is a pre-defined real constant, restr(5, B) returns the restriction of the 
approximate set S (as a function) to B, choose(S', B) returns a subset of B for 
which S is U, and embed(S', S') is an approximate set that is equal to S except 
that it is equal to S' on its domain of definition. 

Up to now we have allowed approximate sets to be arbitrary functions in 
B {T,F, U}. For computer representation we restrict this class to functions 
that are constant on finitely many floating-point boxes. We can represent these 
by a set of boxes for which the approximate set is T, a set of boxes for which 
the approximate set is F, and a set of boxes for which the approximate set is 
U. Now we can easily implement the above functions using this representation. 
Especially, we can implement the computation of f{B) in the refinement function 
by interval methods for computing an overestimation of the range of the function 
/ on a box B ||. 

The only problem lies in the implementation of the propagation function for 
function inversion quantifiers: First, the resulting approximate solution set can 
have an extremely complicated structure, which can be hard or even impossible 
to represent by finitely many floating point boxes. Second, all the approximate 
solution sets occurring in the algorithm have to fulfill the additional property 
of cylindricity [^2|J^ (the clarification of this notion is beyond the scope of this 
paper). To solve these problems remember that the second argument to the 
refinement function is a set B (i.e., a box), whose volume goes to zero during 
the algorithm execution. The smaller this box B is, the more probable it is 
that the argument to the propagation function in the refinement function is an 
approximate solution set that returns the same truth value everywhere. So we 
can remove both problems easily by just returning approximate solution sets that 
return the same truth value everywhere. This means that, given an approximate 
set iS* on a set i3, we define: 

prop[/(^)=^](5) = 

if for aUa;e/(B), S{x) = T then return XxeB.T 
else if for aU x e /(B), S{x) = F then return XxeB.F 
else return \x G B.JJ 

This propagation function always computes constant approximate sets which 
means that it returns the everywhere unknown function very often. In this case 
AQCS will bisect the resulting box and do further calls to the refinement function 



on the single pieces. This need for bisection instead of more inteUigent box 
pruning methods is one of the weaknesses of the current approach. 

5 Example 

We did timings on the following example: 

- f{xi,X2, u, w) = (3a;ia;2 + wu, 2xi + xi) 

- A{x\^X2-,t) —1 < a^i < 1, —1 < a;2 < 1 

- U{u,t) :^ -0.5 < M < 0.5 

- W{w,t) :^ -0.1 < m; < 0.1 

For n = we simply need to solve —1 <x\ < 1,-1 < X2 < 1- For example for 
n = 3, the total number of involved variables is 12. We list timings for n > in 
the table below. All numbers denote seconds needed on a 500MhZ Intel Celeron 
PC running Linux, where oo denotes that the run either needed more than 64MB 
of memory or more than 20min of time. The column with the title "err=0.2" 
shows the time for computing an approximate solution set of error 0.2 on the box 
[—1,1] X [—1,1], the column with the title "single" the time for computing one 
single true box. The last column lists the time for computing an exact symbolic 
solution with the solver QEPCAD |l^. For the latter we did not simply feed 
the input into the program — this would be definitely too hard to solve — but 
we applied the idea of Section ^ also to this case, by computing a solution of 
stage n, and then recursively back-substituting the result for computing earlier 
stages. The figure shows the computed approximate solution set for n = 2 with 
error 0.2 again on the box [—1, 1] x [—1, 1] (green=T, red=F, white=U). 



n err=0.2 single QEPCAD 

1 07 ao lU 

2 19.8 0.3 oo 

3 oo 11.9 oo 

4 00 oo oo 




6 Related Work 



The behavior of various algebraic objects, such as Grobner bases Q or subre- 
sultants [ p^ , under composition is becoming an important topic in computer al- 
gebra. The composition structure of constraints has also been exploited for reuse 
of certain shared interval function evaluation in global optimization p5| , |2l[ | , and 
for computing good range overestimations of functions 113^. 

Inversion of functions on sets is done implicitly by every algorithm for solving 
systems of equations — in this case the input set just contains one zero 
vector. It is mentioned explicitly mostly for computing the solution set of systems 
of inequalities . 

First-order constraints occur frequently in control, and especially robust con- 
trol. Up to now they either have been solved by specialized methods or 
by applying general solvers like QEPCAD In the first case one is usually 
restricted to conditions like linearity, and in the second case one suffers from the 
high run-time complexity of computing exact solutions |4^,0 . We know of only 
one case where general solvers for first-order constraints have been applied to 
discrete-time systems , but they have been frequently applied to continuous 
systems [^,|l5|,0. For non-linear discrete-time systems without perturbations 
or control, interval methods have also proved to be an important tool p^ , |26[ |. 

Apart from the method used in this paper |^ , there there have been several 
successful attempts at solving special cases of first-order constraints, for example 
using classical interval techniques ||3^ , |3^ or constraint satisfaction Q, and very 
often in the context of robust control lll6|,^2|J27|,Ba|. 



7 Conclusion 

We have designed a method for solving composed first-order constraints arising 
in discrete-time robust control by extending a solver (AQCS) developed by one of 
the authors. The result can solve problems that are far too hard to compute for 
state-of-the-art symbolic solvers. However it is still too slow for solving problems 
of practical interest. We believe that the method is promising also in other 
situations where composed first-order constraints occur. 

In problems without function inversion quantifiers, AQCS spends most of 
the time on computing information for atomic constraints, using tightening pof 
and range computation. However, for inputs containing function inversion quan- 
tifiers, this is not the case. Usually more than 90 percent of the time is spent 
on propagating (usually unknown) boxes in tasks like memo lookup or choosing 
boxes for further processing. We envision two basic approaches for dealing with 
this situation: 

— Try to produce less unknown boxes. 

— Try to handle the produced boxes more efficiently. 

For producing less unknown boxes, the following approaches seem promising: 



— Instead of bisection, develop a method similar to tightening [g^l or to box 
consistency methods on the level of function inversion quantifiers, for 
example by implementing the following specification: Given a box B, com- 
pute a box B' such that f{B') C B. Then the fact that B is in the solution 
set of a constraint C{z) implies that B' is in C{f{x)). 

— Pruning boxes that do not contribute to the overall results (in a similar way 
as the monotonicity test in global optimization). 

For handling the produced boxes more efficiently we can: 

— Devise strategies for choosing boxes in the context of function inversion 
quantifiers |33|] . 

— Design efficient algorithms for the memo lookup. 

— Parallelize the method by putting each stage on a different processor. 
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